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Abstract- This paper presents a series of iterative sparse maximum 
likelihood-based approaches (SMLA) with applications to 
synthetic aperture radar (SAR) imaging. By using a particular 
form of Gaussian signal prior, iterative analytical expressions of 
the signal and noise power estimates are obtained by iteratively 
minimizing the stochastic maximum likelihood (SML) function 
with respect to only one scalar parameter at a time, resulting in 
power-based SMLA approaches. However, these power-based 
sparse approaches do not provide the phases of the unknown 
signals. To address this problem, a combined SMLA and 
Maximum A Posteriori (MAP) approach (referred to as the 
SMLA-MAP approach) for estimating the unknown complex- 
valued signals is proposed. The SMLA-MAP derivation is inspired 
by the sparse learning via iterative minimization (SLIM) 
approach, where a modified expression of the SLIM noise power 
estimation is proposed. We also show that SLIM can be viewed as 
a combination of the deterministic ML (DML) and iteratively re- 
weighted least squares (IRLS) approaches. Finally, numerical 
examples of SAR imaging using Slicy data, Backhoe data and 
Gotcha data are generated to compare the performances of the 
proposed and existing approaches . 

Keywords- Synthetic Aperture Radar (SAR); SAR Imaging; 
Sparse Signal Recovery; Maximum Likelihood (ML); Maximum A 
Posteriori (MAP); Sparse Learning via Iterative Minimization 
(SLIM); Iterative Adaptive Approach (IAA). 

I. INTRODUCTION 

Synthetic aperture radar (SAR) has been widely used in a 
variety of applications such as geosciences, remote sensing and 
defense SAR is an active sensor and has all-weather and 
day/night imaging capabilities. There are four distinct modes in 
which a SAR imaging system can operate: scan, stripmap, 
spotlight, and inverse SAR (ISAR). In this paper, we will 
mainly focus on the spotlight mode SAR. In spotlight mode 
SAR, the radar sensor steers its antenna beam to continuously 
illuminate the terrain patch being imaged. It can provide higher 
resolution than the stripmap and scan mode SAR, because it 
views a scene from multiple angles during a single pass [2] . 

There has been considerable interest recently in the use of 
SAR images in automated target recognition and decision- 
making tasks. The success of such tasks depends on how well 
the reconstructed SAR images exhibit certain features of the 
underlying scenes. Among the popular spotlight mode SAR 
image formation algorithms, the polar format algorithm (PFA) 
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, is based on using interpolation and the Fourier transform. 
Another image reconstruction method, suggested by the 
tomographic formulation of SAR [4] , is the filtered 
backprojection (FBP) algorithm [4 ' 5] , where each radial slice in 
the Fourier space is considered as the one dimensional (ID) 
Fourier transform of the projection of the field at the 
corresponding angle. These data-independent algorithms suffer 
from important shortcomings, such as limited resolution, 
speckle, and sidelobe artifacts. Several data-adaptive 
algorithms have been considered for SAR imaging [610] . 
However, these approaches require a large number of data 
snapshots to provide accurate second-order statistics, which is 
hard to satisfy in practice. 

Based on the observation that typical underlying signals of 
interest exhibit a degree of sparsity 1 , Cetin et al. 12] 

proposed a regularization model using £ -norm ( < p < 1 ) 

for enhancement of two dimensional (2D) spotlight mode SAR 
images in the complex-valued image domain. Kragh et al. [13] 

developed another regularized £ p -norm based approach for 

SAR image restoration. These methods, however, like most 
existing sparse signal recovery algorithms, require the delicate 
choice of user parameters. A Bayesian maximum a posteriori 
(MAP) algorithm and a sparse learning via iterative 
minimization (SLIM) algorithm [15] have also been used for 
SAR imaging [21 ' 25] . MAP and SLIM are easy to use and are 
sparse signal recovery algorithms with excellent sidelobe 
suppression and high resolution properties. However, like 
many other sparse signal reconstruction algorithms, MAP and 
SLIM provide biased downward spectral estimates [25] . 
Moreover, although easy to choose, they have a user parameter 
as well. 

Recently, a nonparametric and user parameter-free 
algorithm, referred to as the iterative adaptive approach (IAA) 
[22] , was presented for applications including, for example, SAR 
imaging [20 ' 21 ' 25] and single-antenna radar systems [22] . The 
IAA algorithm provides accurate spectral estimates by 
iteratively updating the IAA co variance matrices, and produces 
SAR images with high resolution and relatively low sidelobes 
[21, 25] 5 mce is computationally expensive, a fast 

implementation of IAA is developed in [23 ' 2 ] to improve its 
computational efficiency by exploiting the Toeplitz structure of 
the IAA co variance matrices. A regularized version of the IAA 
algorithm (IAA-R) is also developed in [20] . However, the 
images generated by IAA are dense while certain recognition 
and decision-making tasks may prefer more sparse images. 



^his means that in an appropriate basis, they can be expressed in terms of a 
small number of nonzero coefficients. 
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In this paper, we first present a series of user parameter-free 
iterative sparse maximum likelihood-based algorithms (SMLA) 
for estimating the sparse signal powers and noise power. We 
then propose a complex-valued SMLA-MAP approach, which 
can not only provide accurate spectral estimation, but also give 
rise to sparse SAR images with high resolution and low 
sidelobe level properties. In Section II, we present briefly the 
complex-valued SAR signal model. In Section III, by using a 
particular form of Gaussian signal prior, iterative analytical 
expressions of the signal and noise power estimates are derived. 
In Section IV, we develop several versions of sparse power- 
based SMLA algorithms. Since the power estimates do not 
contain any information on the phases of the unknown signals, 
in Section V, a combined SMLA and MAP approach (referred 
to as the SMLA-MAP approach) is proposed to estimate the 
unknown complex- valued signals. The SMLA-MAP derivation 
is inspired by the SLIM approach, and a modified expression 
of the SLIM noise power estimation is proposed. We also show 
that the SLIM algorithm can be viewed as a combination of the 
DML and IRLS approaches. In Section VI, numerical examples 
of SAR imaging are generated to compare the performances of 
the proposed and existing approaches, particularly, the well- 
known CoSaMP (Compressive Sampling Matching Pursuit) 
[14]. Finally, conclusions are given in Section VII. 

The following notations are used throughout the paper. 
Matrices and vectors are represented by bold upper and lower 
case characters, respectively. Vectors are by default in column 
orientation, while superscript T , H and * stand for transpose, 
conjugate transpose and complex conjugate, respectively. E(-) , 
Tr(-) and det(-) are the expectation, trace and determinant 
operators, wc(-) is the "vectorization" operator that turns a 
matrix into a vector by stacking the columns of a matrix on top 
of each other, ® denotes the Kronecker matrix product, I is 
the identity matrix, and U m denotes the mth column of I . 

II. PROPLEM FORMULATION 

The phase history data used to form complex-valued SAR 
images can in general be expressed as (e.g., [1 ' 11] ) 



y = Ax+e, 



(1) 



where y e C M is a vector of measurements after appropriate 
pre-processing, e is a complex-valued Gaussian noise vector, 
xgC k is the underlying scene reflectivity (the pixels in a 
scene form a vector), and A e c MxK is a rank deficient 
(M « K) " sensing matrix' ' . 

We assume that the components of the noise vector e are 
independent, zero-mean complex Gaussian random variables 
with unknown variance C . We also assume that e and X are 

def def 

independent, and p = Diag(p) , where p = [p x , . . ., p K ] is the 
signal power vector. The covariance matrix of y that contains 
information about the unknown parameters p and (J is given 
by 

def 

R = APA H + al. 

To form SAR images, we need to estimate the sparse 
vector X from y. The most straightforward method is to 
perform the inverse Fourier transform, which gives the 
conventional delay-and-sum (DAS) or matched filtering 



approach [28] . However, the quality of the DAS image is usually 
unsatisfactory due to the low resolution and high sidelobe 
problems. In this paper, we first propose power -based iterative 
sparse SMLA approaches for estimating p and (J . Owing to 

the complex- valued and potentially random phase nature of the 
reflectivities in SAR images, we then derive a SMLA-MAP 
approach that combines the SMLA and MAP approaches for 
estimating the unknown complex- valued signals X. We note 
that the approaches presented in this paper can also be applied 
to many applications requiring sparse signal recovery 
techniques. 

III. STOCHASTIC ML APPROACH 

In this section, we derive iterative analytical expressions of 
the signal and noise power estimates obtained by iteratively 
minimizing the stochastic maximum likelihood (SML) function 
with respect to only one scalar parameter at a time. Based on 
these expressions, we propose in Section IV the power-based 
SMLA approaches. 

Suppose X and e are circularly Gaussian distributed, and 
y has a circular Gaussian distribution with zero-mean and 
covariance matrix R. Thus, the stochastic negative log- 
likelihood function of y can be expressed as 



L S wl (P><r) = y H R" 1 y + ln(det(R)). 



(2) 



Note that the first term in (2) is a convex function of p and 
(J . However, the second term in (2) can be shown to be a 
concave function of p and (J [ 27 , Theorem 1]. Consequently, 
(2) is a non-convex function of p and <J. Despite this 
difficulty, we prove (see Appendix A) the following result. 

Result 1. The {p k }f =1 and G that minimize (2) can be 

computed iteratively. Assume and d (l) have been 

obtained in the ith iteration, they can be updated at the ( i+1 )th 
iteration as: 



I «H||-l(i) v |2 1 

Pk -(a k H R^a k ) 2+Pk -^^T' k-l,...,K, (3) 



a k H R" i(1) a k 



■ = (y H R" 2(1) y + d (1) Tr(R- 2(1) ) -Tr(R" l(1) ))/Tr(R" 2(1) 



(4) 



where the matrix R can be updated for each iteration as 
R (i) = AP (i) A H +d (i) I and where P (i) = Diag(pf } ,...,p^). 

Through numerical calculations, we have observed that the 

{p k }f =1 and d given by (3) and (4) may be negative; 

therefore, the nonnegativity of the power estimates can be 
enforced at each iteration by setting the negative estimates to 
zero as [ 22 , Eq. (30)], 

p£ +1) =max(0,4 i+1) ) k = l,...,K, 



d (i+1) = max 



(5) 



where cl 1+1> = — - 



la k H R- 1(i) yl 2 



(a k H R- 1(1) a k ) 2 



Fk „H D -l(i) f 



ct = [ y H R y + CT'^T^R )— Tr(R ) ]/Tr(R ). 
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The initialization of p k , k = 1, . . K , can be done with 

the power estimates obtained by means of the delay-and-sum 
(DAS) (or matched filtering) approach, (see, e.g., [28] ) 



A( p) -L^yJl k=i,...,K. 



k,DAS 



la, 



The parameter <T can be initialized for instance as 
d (0) =-llyll 2 . 



(6) 



(7) 



Remark 1. We note that Result 1 can be extended easily to 
the nonuniform white Gaussian noise case where the 
covariance matrix is given by 



def M 



E(e(n)e H (n)) = Diag(a 1? ...,a M ) = Xa m a K+1 X + m 



(8) 



def 



with a K+m = U m , m = 1,. . ,,M, denoting the canonical 
vector with one on its mth entry and zeros elsewhere. Under 
this assumption and from Result 1, the (i+l)th estimate of p k , 
k= 1 , . . . ,K+M, is given by 



- (i+l) _ 1 k 
Pk 



la?R- 1(l) yl 2 



(a^R^a,) 2 a k H R-'% 



1 



where R (i) = AP (i) A H + Y K+M alVa^ and where 



P (1) =Diag(pf\...,p^) 
k = K + l,...,K + M. 



and 



d (i) = D (i) 



IV. SMLA APPROACHES 

As mentioned in Section III, Equations (3) and (4) may 
give negative power estiamtes in the overcomplete basis case 
due to the presence of non-zero terms p k -l/(a"R _1 a k ) and 

a-Tr(R _1 )/Tr(R~ 2 ) . To alleviates this difficulty, we assume 

that 2 Pk =l/(a k I R 1 a k ) and a = Tr(R _1 )/Tr(R" 2 ) , and based 
on Result 1, we propose the following iterative power -based 
SMLA approaches. At the (7+7)th iteration, the signal and 
noise power estimates are computed based on the previous 
iterations. 

A. SMLA-0 Approach 

The estimates of {p k }f =1 and C are updated at the 
(7+7 jth iteration as: 

pP =pf ) (la^R- 1(l) yl 2 ), k = l,...,K, (9) 



v (i+i) _ y K y 



a ' - - 



(10) 



Tr(R" 2(1) ) 



B. SMLA-1 Approach 



The estimates of {p k }f =l and (J are updated at the 
(7+7)th iteration as: 

, _ la?R- 1(i >yl 2 



Pk (a k H R" 1(i >a k ) 2 



Hp-2(i) 

_(i+D _ y K y 



a ' - ■ 



, k = l,...,K, (11) 
(12) 



Tr(R" 2(i) )' 
C. SMLA-2 Approach 

The estimates of {p k }f =1 and <J are updated at the 
(7+7 )th iteration as: 



i)l a t H R- ia V 2 
a«R- ia) a 



Pk "Pk H -l(i) rt ' ~ ' " *' ' KLJ) 



d (i+D = y K y 

Tr(R" 2(i) )' 

D. SMLA- 3 Approach: 

The estimates of {p k }f = i and C are updated at the 
(7+7 jth iteration as: 

(a^R- I(i, a k ) 2 



£(i+D = 



y H R- 2(i) y 

Tr(R- 2<i) )' 



where Rf> = AP, (i) A H +a (i) I 

P^DiagCpf',...,^), 



and 
with 



;<i) 



1 



p-= k = l,...,K. 

In the special case of nonuniform white Gaussian noise 
with covariance matrix given in Remark 2, the SMLA noise 
power estimates must be updated at each iteration as 

q< i+1 >= ' Um I— m = l M (15) 

m (u H R" 1(l) u ) 2 ' ' ' 

V m m / 

where R' 1 ' = AP (i) A H + Y M CT ( ' ) u m u^ and where 

P (i) =Diag( Pl (i) ,...,p«). 

In the following result (proved in Appendix B), we show 
that the SMLA-1 signal and noise power-iterations given by 
(11) and (12) can be obtained by minimizing a weighted least 
squares (WLS) cost function. 

Result 2. The SMLA-1 estimates coincide with the estimates 
p k and d obtained by minimizing the following WLS cost 
function 

def 

p k =argming(p k )=[r y -p k a k ] H C^[r y -p k a k ]. (16) 

Pk 



2 The term p k — l/(a k R a k ) is the estimate of the &-th power 
component given by the Capon estimator [28] . 
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def 



def 



where r y = y <g> y and C' k = C r -pX a k > 

def 

k = l,...,K + l, with C r = R* ®R. 

The implementation steps of these approaches are 
summarized in Table I . 

TABLE I THE SMLA APPROACHES 



Initialization: and (T^ using e.g., (6) and (7) 



Repeat 

. update R (i) = AP (i) A H +a (i) I 
R?> =AP 1 (i) A H +a (i) I. 



,0+1) 



Update p k using SMLA identities (9), (11), (14), or (13). 



• Update <T (?+1) using (12). 



Remark 2. Note that the SMLA-0 power estimate 
expression (9) can be interpreted as 



Pk. 



SMLA-0 



=\x 



k,MAP 



I 2 , k = l,...,K, 



where X kMAP = p k a k R *y is the maximum a posterior 
(MAP) estimate (Eq. (25)) of x k , k = l,... 9 K. 

Comparing SMLA-1 power estimate expression (11) and 
that associated with the IAA approach ([ 22 , Table II]), we see 
that the difference between SMLA-1 and IAA is that the IAA 
power estimate is obtained after the estimation of the 

magnitude of {x k }f =1 (i.e., p kUA =\x kUA I 2 ), by updating 

the IAA covariance matrix APA H for each iteration. 

V. SMLA-MAP 

The main purpose of this section is to derive the SMLA- 
MAP approach based on the power-based SMLA approaches 
presented in Section 4. Our derivation is inspired by the work 
of Wipf et al [16 ' 17] and the derivation of the SLIM approach 

[15] 

A. SLIM 

In this section, we first recall briefly the main iteration 
steps of the SLIM approach [15] . We show that the SLIM 
approach can be viewed as a combination of the deterministic 
ML (DML) and iteratively re-weighted least squares (IRLS) 
approaches. We also quantify the classical SLIM estimate 
expression of the noise power by deriving a new expression 
that depends on the previous estimate of the signal power 
vector. 

The SLIM approach is derived based on the following 
hierarchical Bayesian model [15] : 



where f q (x) is a sparsity-promoting prior for < q < 1 . 

Then, the SLIM estimates of X and (J are obtained by 
minimizing the following MAP criterion: 



def 



{x, a} = arg min L q (x, a) = - log(^(y I x, a)f q (x)f e (a)) 



(19) 



Using (17) and (18), we have 



def 



L q (x,G) = L DML (x,a) + L7(x), 

where L DML (x, a) is the negative log-likelihood function 
for the deterministic ML approach: 

def 1 

L DML (x, a) = Mlog(a) + - II y - Ax II 2 , 

a 

and the regularized function L q e§ (x) is defined as the 
negative of the function f q (x) given by 

def 9 K 

Lj eg (x) =-log(f q (x)) = -£(l x k l q -1). 

q k=i 

Minimizing (19) w.r.t. X is equivalent to minimizing w.r.t. 
X the following function 



Reg. 



y I x, a : N c ( Ax, al), f e (a) oc 1 

def N --(lxJ q -l) 

f q (x) = I> q • 



(17) 



(18) 



1 12 2K 

- II y -Ax II 2 +I* eg (x) = -lly-Axll 2 +-llxll q . 

a a q q 

The solution is given by [ 15 , Eq. (25)] 

x (i+1) =P (i) A H R- 1(i) y, 

with 



(20) 



R (i)d =AP«A H +a (i) I, 

def 

P (i) =Diag(l/p 1 (i) ,...,l/p£ ) ) and pf =\xf r 2 . 

Minimizing (19) w.r.t. <J is equivalent to minimizing the 
DML cost function L DML (.) w.r.t. <J , and the solution is 
given by the following classical SLIM noise power estimate 



1 



£<i+l) ||2 



a^= — lly-Ax 

Therefore, by substituting x <1+1) into (21), we obtain 



(21) 



y-Ax (1+1) =ay, 



where 



=I-(AP W A H )R- ,W =a l0 R" lw . 



def 



(OAHxfi-lCO _^(i)5-Ki) 



Then, the estimate of the noise power can be updated at each 
iteration using the following equation instead of (21): 
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v Hn-2(i) v 

-G+i) _ ^2(i) y K y 



a = a 



M 



(22) 



We note the expression given by (20) is the solution to the 
following £ q optimization problem solved using the iterative 
re- weighted least- squares algorithm (IRLS) (e.g., [18] ) 

^lly-Axll 2 +fX ) lx k | 2 = J^ly-Axll 2 



a 

+ X H P- I(i) x. 



k=l 



The matrix P (l) here can be viewed as the weighting 
matrix. Thus, the SLIM optimization problem (19) is 

equivalent to the following I q optimization problem 



1 



Mlog(cj) + — II y - Ax ll z + II x HJj , 
a 



and its solution is obtained by combing the deterministic ML 
approach to obtain (22) and the IRLS approach to obtain (20). 

B. SM LA-MAP 

We first recall that the SMLA approaches can be viewed as 
approximate solutions of the stochastic ML criterion by 

minimizing the cost function L SML (•) . Note also that the 

estimates of {p k }f =l can be a good tool to determine whether 

=1 are large or small, which is mainly required for 

signal detection applications. However, the powers {p k }f =1 
do not contain any information on the phases of the complex- 
valued signals {x k }f =1 needed by some practical sparse signal 
recovery applications. Moreover, direct calculation of 

{I x k 'lf=i fr° m {Pk }f=i ma y not a ^ wa y s ^ e possible due to 
some unknown scaling factor that is involved (see, e.g., [26] ). To 
alleviate this difficulty, we describe here how to estimate the 
complex-valued signal X and the hyperparameters p and C 

by combining a SMLA approach and the maximum a posteriori 
probability (MAP) approach, resulting in the SMLA-MAP 
method. We note that in the derivation of the SLIM approach, 
the hyperparameters <J and X are estimated by the DML and 
MAP approaches, respectively. Therefore, a question arises as 
well: how can we combine the SML and MAP approaches in 
order to take advantage of the approximate SMLA solution? 

For the special case of Bayesian real- valued data model, 
Wipf et al. have recently proved [ 16 , Theorem 1 and 2] that the 
estimate of X and the hyperparameters (only the signal powers 
are assumed unknown) can be obtained by minimizing an 
auxiliary cost function, constructed by relying on the MAP- 
based cost function (Type I in [16] ) and SML cost function 
(Type II in [16] ). Reminiscent of their approach (e.g., [16 ' 17] ) and 
using the complex-valued data model (1), we propose the 
following auxiliary cost function 



1 



L(x,p,a) = -lly-Axll z +x H P i x + log(det(R)). (23) 
a 



We first assume that p and <J are fixed. Therefore, 

minimizing (23) w.r.t. X is equivalent to minimizing the 
following MAP -based cost function w.r.t. X : 



1 



x = argminL MAP (x) = -lly-Axll / +x n P x. (24) 
x a 



The solution to (24) can be obtained explicitly as 

x M Ap(P^) = PA H R" 1 y- 



(25) 



The next step is to minimize (23) w.r.t. p and <J with 
fixed X. First, inserting the solution x(p,a) into (23) and 
using the following key equality (proved in Appendix C) 

- II y - Ax MAP II 2 +x^ AP P- 1 x MAP = y H R y, (26) 
a 

the cost function (23) becomes the SML cost function 

L,^ (p, a) = y H R" 1 y + log(det(R)). (27) 

As we have shown in Section IV, the SMLA estimates of 
the parameters p and <J given by (9)-(13) and (12) can be 

seen as an approximate solution to the following minimization 
problem 

} = argminL SML (p,a). 

Hence once Psmla anc * ^smla are obtained (at 
convergence of a SMLA approach), and using (25), the final 
estimate of the unknown complex- valued signal X is given by 

def ^_ ^ ^ 

X MAP = X MAP (P SMLA ' °*SMLA ) = ^SMLA^ ^\SMLA y ' 

Where R SMLA = APsMLA A" + ^SMLA 1 !^ > with 

Psmla = Diag(p SMLA ) . 
Consequently, the following result is proved. 
Result 3. Given the final SMLA estimates p SMLA and 

(J SMLA , the MAP estimate of the complex-valued signal X is 
given by 

x map = Psmla A ^Vsmla y* (28) 

From this result, the estimates of the coefficients X and 
the hyperparameters {p, c} may be obtained in a separable 

form. Note also from (28) that if Psmla ^ s K s — sparse 3 , then 

^map ^ s Ks — sparse, with nonzero elements corresponding 

with the 'relevant' basis vectors a k , k = 1,. . K s . 

VI. NUMERICAL AND EXPERIMENTAL EXAMPLES 

In this section, we present the results of the power-based 
SMLA algorithms and the complex-valued SMLA-MAP 
approaches on three SAR data sets: Slicy, Backhoe, and 



3 A vector is said to be if K — sparse if it has at most K nonzero enters. 
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Gotcha. To illustrate the enhanced performance, the proposed 
approaches are compared with the following existing methods: 
FFT, windowed FFT (WFFT), CoSaMP (the user parameter is 
set to be 3000 ) and SLIM (with q = 1 ). Unless otherwise 

specified, the number of iterations 4 is fixed at 10 for all 
iterative algorithms. Note that all images are shown in dB, 
101og 10 for power and 201og 10 for modulus, normalized to 

peak at dB and with a dynamic range of 35 dB. 

A. SAR Images of the Slicy Data 

We use the two-dimensional phase-history data of the Slicy 

target generated at azimuth angle 0° using XPATCH [19] , a 
high frequency electromagnetic scattering prediction code for 
complex 3-D objects. A photo of the Slicy target at 30° 
rotation is shown in Fig. 1. For all obtained images, a 40x40 
center block of the phase-history data is used. 




Fig. 1 Photo of the Slicy object viewed at azimuth angle 45° 

In the first experiment, we show SAR images of the Slicy 
data using the power-based SMLA algorithms presented in 
Section IV. In the second experiment, we show the 
reconstructed Slicy images using the complex-valued SMLA- 
MAP algorithm (i.e., based on the estimates of the complex- 
valued signal) presented in Section V. 

Fig. 2 illustrates the SAR image formed using the WFFT 
approach, where a 2-D Taylor window with a -30 dB sidelobe 
level is applied. We can observe that WFFT has low resolution. 



In contrast, we can see from Fig. 4 that SMLA-0 with 3 
iterations generates a SAR image with much higher resolution 
and lower sidelobes. 



ft 



Fig. 3 2-D SAR image obtained via CoSaMP algorithm 



Fig. 4 2-D SAR image obtained via SMLA-0 algorithm (3 iterations) 

Based on Remark 2, the performance of SMLA-1 and IAA 
are equivalent for high-SNR, where C « . We observe from 
Fig. 5 that the reconstructed SAR image using the SMLA-1 
approach has very high resolution with moderate sidelobe 
levels due to the fact that SMLA-1 is a dense algorithm. 




Fig. 2 2-D SAR image obtained via WFFT algorithm 

Fig. 3 shows the SAR image obtained by CoSaMP, which 
fails to work properly since the object features are deteriorated. 



4 In our simulations, we observe that there is no significant imaging quality 
improvement after 10 iterations. 




Fig. 5 2-D SAR image obtained via SMLA-1 algorithm. 

Figs. 6 and 7 illustrate the imaging performance of the 
SMLA-2 and SMLA-3 approaches, respectively. As can be 
seen in Fig. 6, the SMLA-2 approach attenuates the sidelobe 
further compared with SMLA-1. We conclude that the SMLA- 
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2 approach is sparser than the SMLA-1 approach, and its 
performance is between those of SMLA-0 and SMLA-1. 
Similar to the SMLA-0 approach, the SAR image formed by 
the SMLA-2 approach is obtained with only 3 iterations. It can 
be concluded from Fig. 7 that the SMLA-3 approach produces 
the SAR image with the best appearance among these 
algorithms. However, this is for the noise-free case. 




-35 

Fig. 6 2-D SAR image obtained via SMLA-2 algorithm (3 iterations) 

SMLA-3 



Fig. 7 2-D SAR image obtained via SMLA-3 algorithm 

In the next example, we show in Fig. 8 the modulus of the 
complex-valued SAR image generated using the SLIM-1 
approach (setting q = 1 in (20)). Figs. 9 and 10 demonstrate the 

performance of the proposed complex-valued SMLA-MAP 
approach given by Result 3. Figs. 9 and 10 show the modulus 
of the reconstructed SAR images using the complex-valued 
SMLA-1 -MAP and SMLA-2-MAP approaches, respectively. 
In these figures, the estimates of the complex-valued signal X 
given by SMLA-1 -MAP and SMLA-2-MAP were obtained 
using MAP and the SMLA-1 and SMLA-2 power estimates 
shown in Figs. 5 and 6, respectively. These figures demonstrate 
the effectiveness of these approaches in forming SAR images 
with enhanced resolution and suppressed sidelobes when 
compared to the power-based SMLA-1 and SMLA-2 
approaches. They also show that the SLIM-1 approach is 
somewhat inferior to SMLA-1 -MAP. 



i 



Fig. 8 2-D SAR image obtained via SLEVI-1 algorithm 



Fig. 9 2-D SAR image obtained via SMLA-1 -MAP algorithm 

SMLA-2-MAP 



I 



Fig. 10 2-D SAR image obtained via SMLA-2-MAP algorithm 

In Fig. 1 1 we show the reconstructed image using SMLA-0 
with 10 iterations, where the dominant powers are significantly 
accentuated. This phenomenon is ubiquitous among highly 
sparse algorithms. The reconstructed SMLA-2 image after 10 
iterations is shown in Fig. 12. Fig. 13 shows the number of 
non-zero elements in p , K , versus the iteration number. We 

remark from this figure that after 15 iterations the number of 
non-zero power elements decreases from 1.6 xlO 4 to 91 as the 
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number of iterations increases. Note that both SMLA-0 and 
SMLA-2 are too sparse after 10 iterations. 



To show the robustness of the proposed SMLA-MAP 
algorithms against the presence of noise, we consider further 
the case when the Slicy data is corrupted by white Gaussian 
noise with zero-mean and variance G. Based on the data model 
(1), the SNRis defined as 

def II A x II 2 

SNR = " 11 



Ma 



The SAR images given in Figs. 14-19 are formed with an 
SNR of 5 dB. As can be seen from Figs. 17 and 18, the 
SMLA-1-MAP and SMLA-2-MAP approaches provide better 
SAR images compared with the power-based SMLA-1 and 
SMLA-2 approaches illustrated by Figs. 14 and 15, 
resepectively. We can also see that SMLA-1 -MAP produces a 
better image than SMLA-2-MAP. SMLA-3 and SMLA-3-MAP, 
shown in Figs. 16 and 19, respectively, appear to provide 
poorer images than their SMLA-1 counterparts for this noisy 
data case. 



Fig. 1 1 2-D SAR image (power) obtained via SMLA-0 algorithm after 
convergence (using 10 iterations) 




Fig. 14 2-D SAR image obtained from Slicy data with SNR=5 dB via 
SMLA-1 algorithm 



Fig. 12 2-D SAR image (power) obtained via SMLA-2 algorithm after 
convergence (using 10 iterations) 








15 20 25 

Iteration number 



Fig. 13 The number, K s of non-zero SMLA-0 power estimates versus the 
iteration number 



Fig. 15 2-D SAR image obtained from Slicy data with SNR=5 dB via 
SMLA-2 algorithm (3 iterations) 
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Fig. 16 2-D SAR image obtained from Slicy data with SNR=5 dB via 
SMLA-3 algorithm 



Fig. 17 2-D SAR image obtained from Slicy data with SNR=5 dB via 
SMLA-1-MAP algorithm 



Fig. 19 2-D SAR image obtained from Slicy data with SNR=5 dB via 
SMLA-3 -MAP algorithm 

B. SAR Images of the Backhoe Data 

In this subsection, we present numerical results using the 
Backhoe Data, which is a wideband, full polarization, 
complex-valued backscattered data from a backhoe vehicle in 
free space. Similarly, the Backhoe data is also corrupted by 
white Gaussian noise with zero—mean and variance o . The 3- 
D CAD model of the backhoe and the illustration of the 
backhoe data dome are shown in Figs. 20-21. In our 
experiment, we use the 44x43 HH polarization data with a 
bandwidth of 0.5 GHz centered at 12 GHz. This data is 
collected at a 30° elevation angle and a 5° azimuth cut 
centered around 0° azimuth. 




Fig. 20 The backhoe data: 3-D CAD model 




Fig. 18 2-D SAR image obtained from Slicy data with SNR=5 dB via 
SMLA-2-MAP algorithm (3 iterations) 



kx(GHz) ky(GHz) 

Fig. 21 The backhoe data: K-space data 
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With the same definition of SNR as for the Slicy data, the 
SAR images obtained using Backhoe data are given in Figs. 
22-25 when the SNR is 5 dB. Fig. 22 shows the reconstructed 
SAR image obtained using the FFT approach. We can see from 
this figure that it has poor resolution and high sidelobe level 
problems. Fig. 23 shows the reconstructed image using the 
CoSaMP, which lacks fidelity. Figs. 24-25 illustrate the 
imaging performance of the SMLA-1 and SMLA-1-MAP 
approaches. Again, SMLA-1 -MAP appears to give the best 
performance. 



SMLA-1 -MAP, - AZ: 1deg, ±1.5deg, B: 0.5 GHz 



FFT, - AZ: 1deg, ±1.5deg, B: 0.5 GHz 




Cross Range (meter) 



Fig. 22 2-D SAR image obtained from the Backhoe data set with SNR=5 dB 
via FFT algorithm 



CoSaMP, - AZ: 1deg, ±1.5deg, B: 0.5 GHz 




Cross Range (meter) 



Fig. 23 2-D SAR image obtained from the Backhoe data set with SNR=5 dB 
via CoSaMP algorithm 

SMLA-1, -AZ: 1deg, ±1.5deg, B: 0.5 GHz 









M 











Cross Range (meter) 

Fig. 24 2-D SAR image obtained from the Backhoe data set with SNR=5 dB 
via SMLA-1 algorithm 



• » « / 



Cross Range (meter) 

Fig. 25 2-D SAR image obtained from the Backhoe data set with SNR=5 dB 
via SMLA-1 -MAP algorithm 

C. SAR Images of the Gotcha Data 

We illustrate here the application of SMLA-1 -MAP, which 
has the best performance in the aformentioned numerical 
examples, to circular SAR imaging using the GOTCHA Air 
Force Research laboratory data set. 

The GOTCHA dataset, Volumetric SAR Data Set, Version 
1.0, consists of SAR phase history data collected at X-band 
with a 640 MHz bandwidth with full azimuth coverage at eight 
different elevation angles and full polarization [30] . The imaging 
scene consists of numerous civilian vehicles and calibration 
targets, as shown in Figs. 26-27. Imaging on a small area of the 
illuminated scene may be of more interest in practical 
applications. Below, we aim to reconstruct the image of a 
Chevrolet Malibu in the parking lot from the phase history data 
with full azimuth coverage collected at the first pass for a HH 
polarization channel. 




-50 -40 -30 -20 



10 20 30 40 50 



Fig. 26 2-D SAR image of the GOTCHA scene, from 1 




Fig. 27 Photograph of Chevrolet Malibu from the GOTCHA scene, from 1 
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We use 4° subapertures from 0° to 360° with no overlap, 
which results in a total of 90 subapertures. For each 
subaperture, one 2D spatial image is formed by using a 2D FFT 
operation on the corresponding phase history (k-space) data. 
An 80x80 block of the spatial data centered about the 
Chevrolet Malibu is then chipped out and transformed back 
into k-space using an inverse FFT (IFFT) operation. 

Different spectral analysis algorithms can be applied to the 
so-obtained 80x80 phase history data to get one image for 
each subaperture. With the help of the auxiliary information 
(e.g., the antenna locations, range to scene center, azimuth and 
elevation angles) provided by the GOTCHA data set, each 
image is then projected onto the ground plane and interpolated 
to form a 2D ground image. Finally, all of the 90 2D ground 
images are combined using the noncoherent max magnitude 
operator to yield the reconstructed Malibu image, whose 
dimensions are [5,15] x [-10,0] meters with a grid size of 0.05 
meters in both dimensions. 

Figs. 28-31 shows the reconstructed Chevrolet Malibu SAR 
images obtained by FFT, SMLA-1, SMLA-1-MAP and SLIM- 
1. Again, FFT has low resolution and high sidelobe problems, 
with the front portion of the vehicle mixed in with the 
background, as can be seen from Fig. 28. The Malibu image 
obtained by SMLA-1 is shown in Fig. 29. 

In comparison, we can find that SMLA-1 provides 
significantly higher resolution and lower sidelobe level than 
FFT. The Malibu image obtained via SMLA-1 -MAP is 
illustrated in Fig. 30, from which we can observe that the MAP 
step effectively converts the original dense Malibu images 
formed by SMLA-1 into a much sparser one. The Malibu 
image obtained by SLIM-1 is shown in Fig. 31. 

We observe that the image formed by the SMLA-1 -MAP is 
sparser than that of SLIM-1. It appears that the SMLA-1 -MAP 
satisfactorily balances the tradeoffs between the image 
resolution (note that the SMLA-1 -MAP outperforms SLIM-1 
in terms of resolution) and detail preservation. We conclude 
that the user parameter-free SMLA-1 -MAP again has the best 
performance. 

FFT; HH; Max 



SMLA-1; HH; Max 





Fig. 29 2-D Malibu image obtained via SMLA-1 algorithm 
SMLA-1 -MAP; HH; Max 
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Fig. 30 2-D Malibu image obtained via SMLA-1 -MAP algorithm 
SLIM-1; HH; Max 



S -5 
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Fig. 28 2-D Malibu image obtained via FFT algorithm 



Fig. 3 1 2-D Malibu image obtained via SLIM algorithm 
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VII. CONCLUSIONS 

We have presented a series of power-based iterative SMLA 
approaches for the SAR imaging application based on 
minimizing the SML criterion. Moreover, by exploiting the 
complex-valued nature of the SAR imaging problem, we have 
developed the SMLA-MAP approach, which combines a 
SMLA variant with the MAP approach. Both SMLA and 
SMLA-MAP approaches are user parameter free, making them 
easy to use in practice. Numerical and experimental results on 
various SAR imaging examples demonstrate the merit of the 
SMLA-MAP approach, which provides enhanced resolution 
and suppressed sidelobe levels compared to its SMLA 
counterpart and existing approaches. Among all the competing 
approaches, the newly proposed SMLA- 1 -MAP offers superior 
performance in terms of both SAR imaging quality and 
robustness against the presence of noise. 

APPENDIX A 

PROOF OF RESULT 1 
Define the covariance matrix of the interference and noise 



as 



def 



Q k = R-p k a k a*. 

Applying the matrix inversion lemma, we obtain 

R-^Q^-pAbX. 



(29) 



(30) 



def def 

where b k = Q^a k and p k =(l + p k a k i Q k 1 a k )" 1 . Since 



y H R- 1 y = y H Q- 1 y-p k P k lb k 1 yl 2 , 

ln(det(R)) = ln(det(Q k +p k a k a k H )) 
= ln[(l + p k a k H Q k 1 a k )det(Q k )] 
= ln(det(Q k ))-ln(P k ), 



(31) 



(32) 



where we have used the well-known algebraic identity 
det(I + AB) = det(I + BA). By using a (31) and (32), the 
ML function (2) becomes 

t(p) = ln(det(Q k )) + y "Q-'y - [ln(P k ) + p k P k I b k H y I 2 ] 

Consequently, minimizing (2) w.r.t. p k is equivalent to 
minimizing the following function 



def 



£'( Pk )=m + P k ^Q- k \)-p k P k Ib^yl 2 , (33) 

where ^ k is assumed known and independent of ^ k . 
Note that (33) can be obtained from [22, Eq. (16)] derived 

r = w H 

for a multi-snapshot case by assuming JJ 
Consequently, from [22, Appendix, Eqs. (27) and (28)], the 
unique minimizer of (33) is given by 



Pk = 



a k I Q k 1 (yy H -Q k )Q k " 1 a k 



(a.QkX) 2 



(34) 



Using (30), we have 



a k H Q k a k =Y k (a k 1 R" 1 a k ) 



(35) 



and 



<Qky = y k (a k R"'y), (36) 



def 



where Yk - 1 + Pk a k Q k a k Substituting (35 ) and (36) 
into (34), we obtain the desired expression 

a k R^(yy H -R)R^a k 

P k = -^^kl^V2 ~ + Pk (37) 



(a k R-xr 



1 



la^R-'yl 2 

(afR\) 2+Pr a«R\ 



Differentiating (2) with respect to ® and setting the result 
to zero, we obtain 



Tr(R- 1 (yy H -R)R- 1 ) 



a = 



Tr(R" 2 ) 

and after substituting R — crl for R in the above equation, 
TrCR-^-RJR" 1 ) 



a = 



+ G 



Tr(R 2 ) 

= (y H R" 2 yyrr(R- 2 ) + a - Tr(R _1 )/Tr(R" 2 ). 

Computing Fk and G requires the knowledge of Fk , & , 

and Therefore the algorithm must be implemented 
iteratively as shown in Result 1 . 

APPENDIX B 
PROOF OF RESULT 2 

Differentiating (16) w.r.t. ^ k and setting the result to zero, 
we get (after some simplification) 



(i+ i) = a k^k r y 



(38) 



c 

Applying the inversion lemma to the matrix k , the 
numerator and denominator of Eq. (38) , respectively, can be 
expressed as 



a k H C k " 1 r y =w k (a k H C; 1 r y ), 



and 



a k x a k - w k (a"C/a k ), 



Hr«-1- 



(39) 



(40) 



w, =1+ ' j k 



where 
Thus, 
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(i+ i) a k H C^r y a k H C> y 
Pk o H r^' ] o i- ' K i,...,Jv-hi. ^i; 

a k a k a k a k 

Using the following identities (see e.g., [29, Th. 7.7, 7.16]), 

vec(ABC) = (C T ® A)vec(B), 

( A ® B) ® (C ® D) = AC ® BD, 

the nominator and denominator of Eq. (41) can be expressed, 
respectively, as 

a k H C r _1 r y =\ a *R- l y\ 2 , k = l,...,K, (42) 



and 



a^C; 1 ^ =(a^R- 1 a k ) 2 ^ k = l,...,K, (43) 



and also 



and 



a K H +1 C;V N =y H R- 2 y, 



(44) 



^k+iC'^k+i = Tr(R" 2 ). (45) 

Hence, the result of dividing (42) by (43) gives (11), and 
the result of dividing (44) by (45) gives (12). 

APPENDIX C 
PROOF OF THE KEY EQUALITY 

Since = P A H R ! y , we have 

x1pP _1 Xmap = y H R AP A H R y, (46) 



and 



Then, 



y-Ax MAP =(I-AP A H R 1 )y. (47) 



lly-Ax^ ll 2 =(y-Ax MAP ) H (y-Ax MAP )(48) 
= y H (I - R 1 AP A H )(I - AP A H R 1 )y. 
The left hand side of (26) 

— II y - Ax MAP II 2 +X^ Ap P" 1 X MAP is equal to 



y H [-(I-R AP A H )(I-AP A H R ) + R AP A H R ]y. 

a 

(49) 

If the inside part of the bracket is equal to R , then we 
are done. 



Since 



1 



(I-R AP A H )(I- AP A H R ) 

-Ir 1 



= R [ (R AP A H )(R -AP A H )]R 
a 



= R- 1 [aI]R- 1 , (50) 
the expression (49) becomes 

y H [R (al + AP A H )R ]y = y H [R (R)R ]y = y H R y. 

(51) 
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